# dependencies
library(tidyverse)
library(knitr)
library(kableExtra)
library(boot)
library(parallel)
library(bayestestR)
library(patchwork)
library(mdthemes)
library(lme4)
library(sjPlot)
library(emmeans)
library(ggstance)
library(janitor)
# library(merTools) called via merTools:: to avoid namespace collisions between MASS and dplyr
# set seed for reproducibility
set.seed(42)
# options
options(knitr.table.format = "html") # necessary configuration of tables
# disable scientific notation
options(scipen = 999)
# function to round all numeric vars in a data frame
round_df <- function(df, n_digits = 3) {
require(janitor)
df %>% mutate_if(is.numeric, janitor::round_half_up, digits = n_digits)
}
# create necessary directories
dir.create("../data/processed")
dir.create("../data/results")
#dir.create("models")
# get data
data_trial_level <- read_csv("../data/raw/data_trial_level.csv") %>%
filter(timepoint == "baseline" & (age >= 18 | is.na(age)))
# outliers
data_outliers <- data_trial_level %>%
distinct(unique_id, .keep_all = TRUE) %>%
select(unique_id, domain, mean_rt) %>%
mutate(median_mean_rt = median(mean_rt, na.rm = TRUE),
mad_mean_rt = mad(mean_rt, na.rm = TRUE)) %>%
# exclude median +- 2MAD
mutate(rt_outlier = ifelse(mean_rt < median_mean_rt-mad_mean_rt*2 |
mean_rt > median_mean_rt+mad_mean_rt*2, TRUE, FALSE)) %>%
filter(rt_outlier == FALSE) %>%
select(unique_id, rt_outlier) %>%
full_join(data_trial_level, by = "unique_id") %>%
mutate(rt_outlier = ifelse(is.na(rt_outlier), TRUE, rt_outlier))
data_trimmed <- data_outliers %>%
filter(rt_outlier == FALSE)
# data with confidence intervals
data_estimates_D <- read_csv("../data/processed/data_estimates_D.csv") %>%
filter(method == "bca")
data_estimates_iat_D <- read_csv("../data/processed/data_estimates_iat_D.csv") %>%
mutate(trial_type = "iat",
unique_id = as.factor(unique_id))Most probable estimate among the most probable estimates
data_map_ci_widths <- data_estimates_iat_D %>%
group_by(domain, trial_type) %>%
do(point_estimate(.$ci_width, centrality = "MAP")) %>%
ungroup()
data_map_ci_widths %>%
summarize(map_map = point_estimate(MAP, centrality = "MAP"),
min_map = min(MAP),
max_map = max(MAP)) %>%
unnest(map_map) %>%
rename(MAP_MAP = MAP) %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| MAP_MAP | min_map | max_map |
|---|---|---|
| 0.71 | 0.47 | 0.72 |
By domain and trial type using basic bootstrapping, on the basis that it has the best performance for % of non-zero D scores (further below).
data_map_ci_widths %>%
pivot_wider(names_from = trial_type, values_from = MAP) %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| domain | iat |
|---|---|
| 50 Cent - Britney Spears | 0.70 |
| African Americans - European Americans | 0.70 |
| Artists - Musicians | 0.69 |
| Asians - Whites | 0.71 |
| Astrology - Science | 0.70 |
| Atheism - Religion | 0.71 |
| Athletic People - Intelligent People | 0.72 |
| Avoiding - Approaching | 0.51 |
| Bill Clinton - Hillary Clinton | 0.70 |
| Briefs - Boxers | 0.69 |
| Burger King - McDonald’s | 0.70 |
| Canadian - American | 0.71 |
| Capital Punishment - Imprisonment | 0.71 |
| Career - Family | 0.70 |
| Chaos - Order | 0.56 |
| Coffee - Tea | 0.70 |
| Cold - Hot | 0.71 |
| Conservatives - Liberals | 0.68 |
| Corporations - Nonprofits | 0.71 |
| David Letterman - Jay Leno | 0.71 |
| Denzel Washington - Tom Cruise | 0.69 |
| Determinism - Free will | 0.67 |
| Difficult - Simple | 0.66 |
| Dogs - Cats | 0.66 |
| Dramas - Comedies | 0.69 |
| Drinking - Abstaining | 0.71 |
| Effort - Talent | 0.67 |
| Evolution - Creationism | 0.70 |
| Fat People - Thin People | 0.70 |
| Foreign Places - American Places | 0.71 |
| Friends - Family | 0.71 |
| Gay People - Straight People | 0.71 |
| George Bush - John Kerry | 0.71 |
| Gun Control - Gun Rights | 0.71 |
| Helpers - Leaders | 0.70 |
| Hiphop - Classical | 0.71 |
| Innocence - Wisdom | 0.70 |
| Japan - United States | 0.71 |
| Jazz - Teen Pop | 0.70 |
| Jews - Christians | 0.71 |
| Jocks - Nerds | 0.72 |
| Kobe - Shaq | 0.71 |
| Lawyers - Politicians | 0.71 |
| Lord of the Rings - Harry Potter | 0.70 |
| Manufactured - Natural | 0.63 |
| Meat - Vegetables | 0.69 |
| Meg Ryan - Julia Roberts | 0.70 |
| Microsoft - Apple | 0.71 |
| Money - Love | 0.47 |
| Mother Teresa - Princess Diana | 0.71 |
| Mountains - Ocean | 0.71 |
| Muslims - Jews | 0.70 |
| National Defense - Education | 0.69 |
| New York - California | 0.71 |
| Night - Morning | 0.68 |
| Numbers - Letters | 0.70 |
| Old People - Young People | 0.71 |
| Organized Labor - Management | 0.71 |
| Pants - Skirts | 0.71 |
| Past - Future | 0.65 |
| Pepsi - Coke | 0.71 |
| Poor People - Rich People | 0.56 |
| Private - Public | 0.70 |
| Prolife - Prochoice | 0.71 |
| Protein - Carbohydrates | 0.71 |
| Protestants - Catholics | 0.71 |
| Punishment - Forgiveness | 0.57 |
| Realism - Idealism | 0.71 |
| Reason - Emotions | 0.71 |
| Rebellious - Conforming | 0.68 |
| Receiving - Giving | 0.71 |
| Redsox - Yankees | 0.71 |
| Relaxing - Exercising | 0.71 |
| Republicans - Democrats | 0.68 |
| Rich People - Beautiful People | 0.71 |
| Security - Freedom | 0.71 |
| Single - Married | 0.71 |
| Skeptical - Trusting | 0.51 |
| Solitude - Companionship | 0.70 |
| Southerners - Northerners | 0.70 |
| Speed - Accuracy | 0.71 |
| Stable - Flexible | 0.71 |
| State - Church | 0.71 |
| Strong - Sensitive | 0.71 |
| Tall People - Short People | 0.70 |
| Tax Reductions - Social Programs | 0.71 |
| Team - Individual | 0.71 |
| Technology - Nature | 0.69 |
| Television - Books | 0.70 |
| Tradition - Progress | 0.71 |
| Traditional Values - Feminism | 0.71 |
| Urban - Rural | 0.71 |
| West Coast - East Coast | 0.71 |
| Winter - Summer | 0.71 |
| Wrinkles - Plastic Surgery | 0.71 |
data_ci_width_map_D <- data_estimates_iat_D %>%
group_by(domain, trial_type) %>%
do(point_estimate(.$ci_width, centrality = "MAP")) %>%
ungroup() %>%
mutate(MAP = round_half_up(MAP, 3),
trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4",
trial_type == "iat" ~ "IAT"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4", "IAT")) %>%
mutate(domain = fct_rev(domain))
# # save to disk
# write_csv(data_ci_width_map_D, "../data/results/data_ci_width_map_D.csv")
# plot
p_ci_widths <-
ggplot(data_ci_width_map_D, aes(MAP, domain)) +
geom_point(position = position_dodge(width = 0.8)) +
mdthemes::md_theme_linedraw() +
facet_wrap(~ trial_type, ncol = 4, nrow = 1) +
labs(x = "Highest probability (MAP) 95% CI width",
y = "") +
theme(legend.position = "top")
p_ci_widthsp_cis_by_domain <-
data_estimates_iat_D %>%
arrange(estimate) %>%
group_by(domain) %>%
mutate(ordered_id = row_number()/n()) %>%
ungroup() %>%
ggplot() +
geom_linerange(aes(x = ordered_id, ymin = ci_lower, ymax = ci_upper, color = sig),
alpha = 1) +
geom_point(aes(ordered_id, estimate), size = 0.5, shape = "square") +
geom_hline(yintercept = 0, linetype = "dotted") +
mdthemes::md_theme_linedraw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "top") +
scale_color_viridis_d(end = 0.6, direction = -1) +
xlab("Ranked participant") +
ylab("*D* score") +
labs(color = "95% CI excludes zero point") +
facet_wrap(~ domain, ncol = 6)
p_cis_by_domaindata_diff_zero <-
bind_rows(
mutate(data_estimates_D, measure = "IRAP"),
mutate(data_estimates_iat_D, measure = "IAT")
) %>%
mutate(measure = fct_relevel(measure, "IRAP", "IAT"),
domain = as.factor(domain),
trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4",
trial_type == "iat" ~ "IAT"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4", "IAT")) %>%
group_by(domain, trial_type, measure) %>%
summarize(proportion_diff_zero = mean(sig),
variance = plotrix::std.error(sig)^2,
.groups = "drop") %>%
# model cannot be run on zero variance or 0 or 1 logit, so offset by a minuscule amount
mutate(proportion_diff_zero_temp = case_when(proportion_diff_zero < 0.001 ~ 0.001,
proportion_diff_zero > 0.999 ~ 0.999,
TRUE ~ proportion_diff_zero),
proportion_diff_zero_logit = boot::logit(proportion_diff_zero_temp)) %>%
select(-proportion_diff_zero_temp) %>%
filter(!(proportion_diff_zero == 0 & variance == 0)) %>%
mutate(variance = ifelse(variance == 0, 0.0001, variance))
data_diff_zero## # A tibble: 227 × 6
## domain trial…¹ measure propo…² varia…³ propo…⁴
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 50 Cent - Britney Spears IAT IAT 0.46 2.51e-3 -0.160
## 2 African Americans - European America… IAT IAT 0.46 2.51e-3 -0.160
## 3 Artists - Musicians IAT IAT 0.41 2.44e-3 -0.364
## 4 Asians - Whites IAT IAT 0.48 2.52e-3 -0.0800
## 5 Astrology - Science IAT IAT 0.5 2.53e-3 0
## 6 Atheism - Religion IAT IAT 0.58 2.46e-3 0.323
## 7 Athletic People - Intelligent People IAT IAT 0.49 2.52e-3 -0.0400
## 8 Avoiding - Approaching IAT IAT 0.95 4.80e-4 2.94
## 9 Bill Clinton - Hillary Clinton IAT IAT 0.43 2.48e-3 -0.282
## 10 Body image Trial … IRAP 0.190 7.71e-3 -1.45
## # … with 217 more rows, and abbreviated variable names ¹trial_type,
## # ²proportion_diff_zero, ³variance, ⁴proportion_diff_zero_logit
## # ℹ Use `print(n = ...)` to see more rows
# # save to disk
# write_csv(data_diff_zero, "../data/results/data_diff_zero.csv")p_diff_zero <-
data_diff_zero %>%
filter(measure == "IAT") %>%
mutate(domain = fct_rev(factor(domain))) %>%
ggplot(aes(proportion_diff_zero, domain)) +
geom_linerangeh(aes(xmin = proportion_diff_zero - sqrt(variance)*1.96,
xmax = proportion_diff_zero + sqrt(variance)*1.96),
position = position_dodge(width = 0.75)) +
geom_point(position = position_dodge(width = 0.75)) +
#scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
#scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
mdthemes::md_theme_linedraw() +
labs(x = "Proportion of scores different from zero point",
y = "") +
theme(legend.position = "top")
p_diff_zero# fit model
fit_diff_zero <-
lmer(proportion_diff_zero_logit ~ 1 + measure + (1 | domain/trial_type),
weights = 1/variance,
data = data_diff_zero,
# solution from https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
control = lmerControl(check.nobs.vs.nlev = "ignore",
check.nobs.vs.nRE = "ignore"))
# extract marginal means
results_emm_diff_zero <-
summary(emmeans(fit_diff_zero, ~ measure)) %>%
dplyr::select(measure, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL)
# extract re Tau
results_re_tau_diff_zero <- fit_diff_zero %>%
merTools::REsdExtract() %>%
as_tibble(rownames = "re") %>%
rename(tau = value)
# combine
results_diff_zero <- results_emm_diff_zero %>%
mutate(pi_lower = estimate - (1.96 * sqrt(se^2 + results_re_tau_diff_zero$tau[1]^2)), # as in metafor package's implementation of prediction intervals, see metafor::predict.rma.R
pi_upper = estimate + (1.96 * sqrt(se^2 + results_re_tau_diff_zero$tau[1]^2))) |>
select(-se) |>
mutate_if(is.numeric, boot::inv.logit)
# plot
p_prop_nonzero <-
ggplot(results_diff_zero, aes(measure, estimate)) +
geom_linerange(aes(ymin = pi_lower, ymax = pi_upper), size = 0.5, position = position_dodge(width = 0.8), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1.25, position = position_dodge(width = 0.8)) +
geom_point(position = position_dodge(width = 0.8), size = 2.5) +
mdthemes::md_theme_linedraw() +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
#scale_color_viridis_d(alpha = 1, begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
#scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
scale_x_discrete(labels = c("IRAP D scores", "IAT D scores")) +
labs(x = "",
y = "Proportion of participants with non-zero scores<br/>") +
theme(legend.position = "none") +
coord_flip(ylim = c(0, 1))
p_prop_nonzeroresults_diff_zero %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| measure | estimate | ci_lower | ci_upper | pi_lower | pi_upper |
|---|---|---|---|---|---|
| IRAP | 0.14 | 0.12 | 0.16 | 0.04 | 0.39 |
| IAT | 0.54 | 0.51 | 0.58 | 0.23 | 0.83 |
# tests
data_emms_diff_zero <- emmeans(fit_diff_zero, list(pairwise ~ measure), adjust = "holm")
summary(data_emms_diff_zero)$`pairwise differences of measure` %>%
as.data.frame() %>%
select(comparison = 1, p.value) %>%
mutate(p.value = ifelse(p.value < .001, "< .001", round_half_up(p.value, 3))) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| comparison | p.value |
|---|---|
| IRAP - IAT | < .001 |
Within domain and trial type.
Note: Discriminability between a score and zero can be determined using the CI, because zero is a known value and only the score is measured with uncertainty. However, discriminability between two scores must take into account the uncertainty in the estimation of both scores. Weir (2005) argues that such an interval can be estimated by expanding the CIs by sqrt(2). Here I refer to these intervals as Discriminability Intervals (DIs).
# helper function to apply workflow to each resample
discriminability <- function(data, i) {
data_with_indexes <- data[i,] # boot function requires data and index
estimate <- data_with_indexes$estimate
di_lower <- data_with_indexes$di_lower
di_upper <- data_with_indexes$di_upper
n_estimate <- length(estimate)
n_di_lower <- length(di_lower)
n_di_upper <- length(di_upper)
r_estimate <- sum(rank(c(estimate, di_lower))[1:n_estimate])
r_di_upper <- sum(rank(c(di_upper, estimate))[1:n_di_upper])
prob_estimate_inferior_to_di_lower <- 1 - (r_estimate / n_estimate - (n_estimate + 1) / 2) / n_di_lower
prob_estimate_superior_to_di_upper <- 1 - (r_di_upper / n_di_upper - (n_di_upper + 1) / 2) / n_estimate
probability_estimates_outside_cis <- (prob_estimate_inferior_to_di_lower + prob_estimate_superior_to_di_upper)
return(probability_estimates_outside_cis)
}
bootstrap_discriminability <- function(data){
require(dplyr)
require(boot)
fit <-
boot::boot(data = data,
statistic = discriminability,
R = 5000,
sim = "ordinary",
stype = "i",
parallel = "multicore",
ncpus = parallel::detectCores()-1)
results <- boot::boot.ci(fit, conf = 0.95, type = "bca")
output <-
tibble(
estimate = fit$t0,
ci_lower = results$bca[4],
ci_upper = results$bca[5]
)
return(output)
}
# irap data
data_discriminability_D <- read_csv("../data/results/data_discriminability_D.csv") %>%
filter(method == "bca")
# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("../data/results/data_discriminability_iat_D.csv")) {
data_discriminability_iat_D <- read_csv("../data/results/data_discriminability_iat_D.csv")
} else {
# bootstrap D scores
data_discriminability_iat_D <- data_estimates_iat_D |>
# expand CIs by sqrt(2) to form discriminability intervals
mutate(di_lower = estimate - (estimate - ci_lower)*sqrt(2),
di_upper = estimate + (ci_upper - estimate)*sqrt(2)) |>
select(unique_id, domain, trial_type, estimate, di_upper, di_lower) |>
group_by(domain, trial_type) |>
do(bootstrap_discriminability(data = .)) |>
ungroup() |>
rename(proportion_discriminable = estimate) |>
mutate(variance = ((ci_upper - ci_lower)/(1.96*2))^2,
domain = as.factor(domain),
trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4", "iat"),
measure = "IAT")
# save to disk
write_csv(data_discriminability_iat_D, "../data/results/data_discriminability_iat_D.csv")
}# combine
data_discriminability_combined <-
bind_rows(
mutate(data_discriminability_D, measure = "IRAP"),
mutate(data_discriminability_iat_D, measure = "IAT")
) %>%
mutate(measure = fct_relevel(measure, "IRAP", "IAT"),
trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4",
trial_type == "iat" ~ "IAT"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4", "IAT")) %>%
filter(!(proportion_discriminable == 0 & variance == 0)) %>%
mutate(variance = ifelse(variance == 0, 0.0001, variance)) |>
# model cannot be run on zero variance or 0 or 1 logit, so offset by a minuscule amount
mutate(
proportion_discriminable_temp = case_when(proportion_discriminable < 0.001 ~ 0.001,
proportion_discriminable > 0.999 ~ 0.999,
TRUE ~ proportion_discriminable),
proportion_discriminable_logit = boot::logit(proportion_discriminable_temp)
) %>%
select(-proportion_discriminable_temp)
p_discriminability <-
data_discriminability_combined %>%
filter(measure == "IAT") %>%
mutate(domain = fct_rev(factor(domain))) %>%
ggplot(aes(proportion_discriminable, domain)) +
geom_linerangeh(aes(xmin = proportion_discriminable - sqrt(variance)*1.96,
xmax = proportion_discriminable + sqrt(variance)*1.96),
position = position_dodge(width = 0.75)) +
geom_point(position = position_dodge(width = 0.75)) +
#scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
#scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
mdthemes::md_theme_linedraw() +
facet_wrap(~ trial_type, ncol = 4) +
labs(x = "Proportion of participants<br/>whose scores differ from one another<br/>",
y = "",
color = "Scoring method",
shape = "Scoring method") +
theme(legend.position = "top")
p_discriminability# fit meta analytic model
fit_disciminability <-
lmer(proportion_discriminable_logit ~ 1 + measure + (1 | domain/trial_type),
weights = 1/variance,
data = data_discriminability_combined,
# solution from https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
control = lmerControl(check.nobs.vs.nlev = "ignore",
check.nobs.vs.nRE = "ignore"))
# extract marginal means
results_emm_disciminability <-
summary(emmeans(fit_disciminability, ~ measure)) %>%
dplyr::select(measure, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL)
# extract re Tau
results_re_tau_disciminability <- fit_disciminability %>%
merTools::REsdExtract() %>%
as_tibble(rownames = "re") %>%
rename(tau = value)
# combine
results_disciminability <- results_emm_disciminability %>%
mutate(pi_lower = estimate - (1.96 * sqrt(se^2 + results_re_tau_disciminability$tau[1]^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
pi_upper = estimate + (1.96 * sqrt(se^2 + results_re_tau_disciminability$tau[1]^2))) |>
select(-se) |>
mutate_if(is.numeric, boot::inv.logit)
# plot
p_prop_discriminable <-
ggplot(results_disciminability, aes(measure, estimate)) +
geom_linerange(aes(ymin = pi_lower, ymax = pi_upper), size = 0.5, position = position_dodge(width = 0.8), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1.25, position = position_dodge(width = 0.8)) +
geom_point(position = position_dodge(width = 0.8), size = 2.5) +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
#scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
#scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
scale_x_discrete(labels = c("IRAP D scores", "IAT D scores")) +
mdthemes::md_theme_linedraw() +
labs(x = "",
y = "Proportion of participants<br/>whose scores differ from one another<br/>") +
theme(legend.position = "none") +
coord_flip(ylim = c(0, 1))
p_prop_discriminable results_disciminability %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| measure | estimate | ci_lower | ci_upper | pi_lower | pi_upper |
|---|---|---|---|---|---|
| IRAP | 0.07 | 0.06 | 0.08 | 0.03 | 0.17 |
| IAT | 0.43 | 0.40 | 0.46 | 0.21 | 0.67 |
# tests
data_emms_disciminability <- emmeans(fit_disciminability, list(pairwise ~ measure), adjust = "holm")
summary(data_emms_disciminability)$`pairwise differences of measure` %>%
as.data.frame() %>%
select(comparison = 1, p.value) %>%
mutate(p.value = ifelse(p.value < .001, "< .001", round_half_up(p.value, 3))) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| comparison | p.value |
|---|---|
| IRAP - IAT | < .001 |
NB observed range of confidence intervals
## calculate observed ranges
observed_range_estimates_D <- data_estimates_D %>%
group_by(domain, trial_type) %>%
dplyr::summarize(min = min(ci_lower, na.rm = TRUE),
max = max(ci_upper, na.rm = TRUE),
.groups = "drop") %>%
mutate(range = max - min)
observed_range_estimates_iat_D <- data_estimates_iat_D %>%
group_by(domain) %>%
dplyr::summarize(min = min(ci_lower, na.rm = TRUE),
max = max(ci_upper, na.rm = TRUE),
.groups = "drop") %>%
mutate(range = max - min)
# calculate CI / range
data_ci_width_proportions_D <- data_estimates_D %>%
# join this data into the original data
full_join(observed_range_estimates_D, by = c("domain", "trial_type")) %>%
# calculate ci width as a proportion of observed range
mutate(ci_width_proportion = ci_width / range) %>%
mutate(measure = "IRAP")
data_ci_width_proportions_iat_D <- data_estimates_iat_D %>%
# join this data into the original data
full_join(observed_range_estimates_iat_D, by = "domain") %>%
# calculate ci width as a proportion of observed range
mutate(ci_width_proportion = ci_width / range) %>%
mutate(measure = "IAT")
# combine
data_ci_width_proportions_combined <-
bind_rows(
data_ci_width_proportions_D,
data_ci_width_proportions_iat_D
) %>%
mutate(measure = fct_relevel(measure, "IRAP", "IAT"),
domain = as.factor(domain),
trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4", "iat")) %>%
# logit transform
mutate(ci_width_proportion_temp = case_when(ci_width_proportion < 0.0001 ~ 0.0001,
ci_width_proportion > 0.9999 ~ 0.9999,
TRUE ~ ci_width_proportion),
ci_width_proportion_logit = boot::logit(ci_width_proportion_temp)) %>%
select(-ci_width_proportion_temp)# fit model
fit_ci_width_proportions <-
lmer(ci_width_proportion_logit ~ 1 + measure + (1 | domain/trial_type),
data = data_ci_width_proportions_combined,
# solution from https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
control = lmerControl(check.nobs.vs.nlev = "ignore",
check.nobs.vs.nRE = "ignore"))
# extract marginal means
results_emm_ci_width_proportions <-
summary(emmeans(fit_ci_width_proportions, ~ measure)) %>%
dplyr::select(measure, estimate = emmean, se = SE, ci_lower = asymp.LCL, ci_upper = asymp.UCL)
# extract re Tau
results_re_tau_ci_width_proportions <-
merTools::REsdExtract(fit_ci_width_proportions) %>%
as_tibble(rownames = "re") %>%
rename(tau = value)
# combine
results_ci_width_proportions <- results_emm_ci_width_proportions %>%
mutate(pi_lower = estimate - (1.96 * sqrt(se^2 + results_re_tau_ci_width_proportions$tau[1]^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
pi_upper = estimate + (1.96 * sqrt(se^2 + results_re_tau_ci_width_proportions$tau[1]^2))) %>%
select(-se) %>%
mutate_if(is.numeric, boot::inv.logit)
# plot
p_ci_width_proportion_observed_range <-
ggplot(results_ci_width_proportions, aes(measure, estimate,
)) +
geom_linerange(aes(ymin = pi_lower, ymax = pi_upper), size = 0.5, position = position_dodge(width = 0.8), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1.25, position = position_dodge(width = 0.8)) +
geom_point(position = position_dodge(width = 0.8), size = 2.5) +
#scale_shape_discrete(labels = c("IRAP", "IAT")) +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Better)", "0.25", "0.50", "0.75", "1.00<br/>(Worse)")) +
#scale_shape_manual(labels = c("IRAP", "IAT"), values = c(15, 16)) +
#scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("IRAP", "IAT")) +
scale_x_discrete(labels = c("IRAP D scores", "IAT D scores")) +
mdthemes::md_theme_linedraw() +
labs(x = "",
y = "Proportion of observed range covered<br/>by individual participants' 95% CIs") +
theme(legend.position = "none") +
coord_flip(ylim = c(0, 1))
p_ci_width_proportion_observed_rangeresults_ci_width_proportions %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| measure | estimate | ci_lower | ci_upper | pi_lower | pi_upper |
|---|---|---|---|---|---|
| IRAP | 0.52 | 0.50 | 0.53 | 0.43 | 0.60 |
| IAT | 0.26 | 0.26 | 0.27 | 0.21 | 0.33 |
# tests
data_emms_ci_width_proportions <- emmeans(fit_ci_width_proportions, list(pairwise ~ measure), adjust = "holm")
summary(data_emms_ci_width_proportions)$`pairwise differences of measure` %>%
as.data.frame() %>%
select(comparison = 1, p.value) %>%
mutate(p.value = ifelse(p.value < .001, "< .001", round_half_up(p.value, 3))) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| comparison | p.value |
|---|---|
| IRAP - IAT | < .001 |
p_cis_by_domainggsave(filename = "plots/figure_4_cis_by_domain_irap_vs_iat.pdf",
plot = p_cis_by_domain,
device = "pdf",
# path = NULL,
# dpi = 300,
units = "in",
width = 8,
height = 16,
limitsize = TRUE)p_ci_widthsggsave(filename = "plots/figure_5_ci_widths_irap_vs_iat.pdf",
plot = p_ci_widths,
device = "pdf",
# path = NULL,
# dpi = 300,
units = "in",
width = 6,
height = 12,
limitsize = TRUE)p_combined <-
p_prop_nonzero +
p_prop_discriminable +
p_ci_width_proportion_observed_range +
plot_layout(ncol = 1) #, guides = "collect") & theme(legend.position = "top")
p_combinedggsave(filename = "plots/figure_6_metaanalyses_irap_vs_iat.pdf",
plot = p_combined,
device = "pdf",
# path = NULL,
# dpi = 300,
units = "in",
width = 5,
height = 5,
limitsize = TRUE)sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_IE.UTF-8/en_IE.UTF-8/en_IE.UTF-8/C/en_IE.UTF-8/en_IE.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] janitor_2.1.0 ggstance_0.3.5 emmeans_1.7.3 sjPlot_2.8.10
## [5] lme4_1.1-29 Matrix_1.4-1 mdthemes_0.1.0 patchwork_1.1.1
## [9] bayestestR_0.12.1 boot_1.3-28 kableExtra_1.3.4 knitr_1.39
## [13] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4
## [17] readr_2.1.2 tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6
## [21] tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.1-1 minqa_1.2.4 colorspace_2.0-3
## [4] ellipsis_0.3.2 sjlabelled_1.2.0 estimability_1.3
## [7] snakecase_0.11.0 markdown_1.1 parameters_0.18.1
## [10] fs_1.5.2 gridtext_0.1.4 ggtext_0.1.1
## [13] rstudioapi_0.13 listenv_0.8.0 furrr_0.3.0
## [16] farver_2.1.1 bit64_4.0.5 fansi_1.0.3
## [19] mvtnorm_1.1-3 lubridate_1.8.0 xml2_1.3.3
## [22] codetools_0.2-18 splines_4.2.0 sjmisc_2.8.9
## [25] jsonlite_1.8.0 nloptr_2.0.3 ggeffects_1.1.2
## [28] pbkrtest_0.5.1 broom_0.8.0 dbplyr_2.1.1
## [31] broom.mixed_0.2.9.4 shiny_1.7.1 effectsize_0.6.0.1
## [34] compiler_4.2.0 httr_1.4.3 sjstats_0.18.1
## [37] backports_1.4.1 assertthat_0.2.1 fastmap_1.1.0
## [40] cli_3.3.0 later_1.3.0 htmltools_0.5.2
## [43] tools_4.2.0 coda_0.19-4 gtable_0.3.0
## [46] glue_1.6.2 merTools_0.5.2 Rcpp_1.0.9
## [49] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.4.1
## [52] svglite_2.1.0 nlme_3.1-157 iterators_1.0.14
## [55] insight_0.18.0 xfun_0.31 globals_0.14.0
## [58] rvest_1.0.2 mime_0.12 lifecycle_1.0.1
## [61] future_1.25.0 MASS_7.3-56 zoo_1.8-10
## [64] scales_1.2.0 vroom_1.5.7 promises_1.2.0.1
## [67] hms_1.1.1 sandwich_3.0-1 yaml_2.3.5
## [70] sass_0.4.1 stringi_1.7.8 highr_0.9
## [73] foreach_1.5.2 plotrix_3.8-2 blme_1.0-5
## [76] rlang_1.0.4 pkgconfig_2.0.3 systemfonts_1.0.4
## [79] arm_1.12-2 evaluate_0.15 lattice_0.20-45
## [82] labeling_0.4.2 bit_4.0.4 tidyselect_1.1.2
## [85] parallelly_1.31.1 magrittr_2.0.3 R6_2.5.1
## [88] generics_0.1.2 multcomp_1.4-19 DBI_1.1.2
## [91] pillar_1.8.0 haven_2.5.0 withr_2.5.0
## [94] abind_1.4-5 survival_3.3-1 datawizard_0.4.1
## [97] performance_0.9.1 modelr_0.1.8 crayon_1.5.1
## [100] utf8_1.2.2 tzdb_0.3.0 rmarkdown_2.14
## [103] grid_4.2.0 readxl_1.4.0 reprex_2.0.1
## [106] digest_0.6.29 webshot_0.5.3 xtable_1.8-4
## [109] httpuv_1.6.5 munsell_0.5.0 viridisLite_0.4.0
## [112] bslib_0.3.1